Method and apparatus for segmenting small structures in images

ABSTRACT

A method for segmenting a small feature in a multidimensional digital array of intensity values in a data processor computes an edge metric along each ray of a plurality of multidimensional rays originating at a local intensity extreme (local maximum or minimum). A multidimensional point corresponding to a maximum edge metric on each said ray is identified as a ray edge point. Every point on each ray from the local extreme to the ray edge point is labeled as part of the small object. Further points on the feature are grown by labeling an unlabeled point if the unlabeled point is adjacent to a labeled point, and the unlabeled point has a more extreme intensity than the labeled point, and the unlabeled point is closer than the labeled point to the local extreme. The resulting segmentation is quick, and identifies boundaries of small features analogous to boundaries identified by human analysts, and does not require statistical parameterizations or thresholds manually determined by a user.

Notice: More than one reissue application has been filed for the reissueof U.S. Pat. No. 7,106,893. The reissue applications are applicationSer. No. 12/210,107, which was filed on Sep. 12, 2008 (the presentapplication), and application Ser. No. 13/314,021, which was filed onDec. 7, 2011 and is a continuation of the present application.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a broadening reissue of U.S. Pat. No. 7,106,893,issued Sep. 12, 2006, from U.S. patent application Ser. No. 10/716,797,filed Nov. 18, 2003, which is a continuation of U.S. patent applicationSer. No. 09/305,018, filed May 4, 1999, now abandoned, which claims thebenefit of provisional patent application Ser. No. 60/084,125 filed onMay 4, 1998, the entire disclosure of which is incorporated herein byreference all of which are incorporated herein by reference in theirentireties.

FIELD OF THE INVENTION

The present invention relates to data processing of intensity dataarranged in a multidimensional array. More particularly, the inventionrelates to a method, an apparatus, and computer program products forrapidly segmenting multidimensional intensity data by which points inone or more small structures contained in the data are labeled.

BACKGROUND OF THE INVENTION

Digital imagery and other multidimensional digital arrays of intensityare routinely collected using digital sensors and arrays of chargecoupled devices (CCDs). The resulting data arrays are analyzed todetermine patterns and detect features in the data. For example, colorimages of a battle scene are analyzed to detect targets, and radiographsand sonograms of human and animal bodies are examined to detect tumorsand other indications of injury or disease. As the number and complexityof these digital data arrays to be analyzed increase or the timerequired to perform the analyses decreases, automated and machineassisted analysis becomes more critical. Some statistically basedautomated procedures for detecting features in a multidimensional arrayare adequate when the feature encompasses many points in the array, i.e.when the feature is large, but fail to perform well as the feature to bedetected becomes small. Some procedures perform well when tuned to aparticular problem through experimental adjustment of many parameters,but such tuning may place an undue burden on time limited or experiencelimited personnel. Typical problems encountered with such automatedanalysis of small structures in multidimensional arrays are illustratedfor the case of automatic detection of microcalcification candidates inmammograms.

Breast cancer has the highest incidence among all cancer types inAmerican women, causing 1 woman in 8 to develop the disease in herlifetime. Every year, about 182,000 new cases are diagnosed with breastcancer and about 46,000 women die of this disease. The 5-year survivalfor women with breast cancer improves significantly with early diagnosisand treatment. To enable early detection, the American Cancer Society(ACS) recommends a baseline mammogram for all women by the age of 40, amammogram approximately every other year between the ages of 40 and 50,and a mammogram every year after the age of 50. It is possible that thevolume of mammography will become one of the highest among clinicalX-ray procedures since more than 30 million women in the U.S. are abovethe age of 50 and 41% are known to follow the ACS guidelines.

Besides the volume problem, an additional difficulty of early detectionof breast cancer in mammograms is the subtlety of the early signal. Amicrocalcification cluster, an early sign of breast cancer that maywarrant biopsy, is commonly defined as three or more microcalcificationspresent in 1 cm² on a mammogram. These clusters are often difficult todetect due to their small size and their similarity to other tissuestructures. The width of an individual microcalcification is less than 2mm. The etiology of microcalcifications includes lobular, ductal orepithelial hyperplasia, secretion of calcium salts by epithelial cells,adenosis, as well as calcification of necrotic debris due to carcinoma.Up to 50% of breast cancer cases exhibit microcalcification clusters,and 20-35% of clusters in the absence of a mass are related to malignantgrowth. In many cases a cluster is the first and only sign that allowstimely intervention.

The increasing pressure to interpret large numbers of mammograms and thesubtlety of many early signs increase the likelihood of missing breastcancer. A reliable automated system that indicates suspicious structuresin mammograms can allow the radiologist to focus rapidly on the relevantparts of the mammogram and it can increase the effectiveness andefficiency of radiology clinics. In the detection of breast cancer,false negatives may cause a delay in the diagnosis and treatment of thedisease while false positives cause unwarranted biopsy examinations.Therefore, both sensitivity and specificity need to be maximized, with arelatively higher priority on sensitivity, which has a more vital role.

A common approach used for detecting microcalcifications in mammogramsstarts by segmenting candidate structures and subsequently applyingfeature extraction and pattern recognition to distinguishmicrocalcifications from background tissue among the candidates. In thisprocess, segmentation plays an essential role since the quantitativefeatures that represent each candidate structure, such as size,contrast, and sharpness, depend on the region indicated by segmentation.Furthermore, to process all possible candidate structures, aconsiderably large number of background structures need to be segmented,making fast segmentation desirable.

Several techniques for segmentation have been applied tomicrocalcifications. One segmentation technique is based on localthresholding for individual pixels using the mean pixel value and rootmean square (rms) noise fluctuation in a selected region around thethresholded pixel. The threshold for a pixel is set as the mean valueplus the rms noise value multiplied by a selected coefficient. Astructure is segmented by connecting pixels that exceed the threshold.Both parameters that have to be selected, size of region and thresholdcoefficient, are critical to this method. If a microcalcification isclose to another microcalcification or bright structure, the window usedto compute the rms noise value around the first microcalcification willinclude the other bright structures, and the noise rms may beoverestimated, thus setting the threshold too high. On the other hand,if the selected region is too small, it will not contain sufficientbackground pixels when placed on large microcalcifications.

Such a window size needs to be selected in a second segmentationalgorithm as well, where local thresholding is used by setting athreshold for small square sub images. The threshold is based on anexpected bimodal intensity distribution in a window of selected sizethat contains the sub-image to be segmented. If the distribution is notbimodal, then the threshold is set by using 5 different positions of thewindow each containing the sub-image to be segmented. The existence of abimodal distribution in at least one window is essential for thisalgorithm.

Other segmentation methods start with seed pixels and grow a region byadding pixels. They also require selection of a window size andthreshold parameters. The localized implementation of region growingdepends on the selected window size and the threshold for absolutedifference in gray level between the seed pixel and a pixel to be addedto the region.

One segmentation algorithm uses several steps that include high-passfiltering, difference of Gaussinan filtering, four computations of thestandard deviation of the image, a smoothing, an opening, as well as aniterative thickening process with two erosions, two intersections and aunion operation in each iteration. More than ten parameters have to beselected, including widths of Gaussian distributions, thresholdcoefficients, and diameters of morphological filtering elements.

A segmentation algorithm that operates without parametric distributionmodels, local statistics windows, or manually adjustable thresholds isdesirable.

A segmentation method that is fast is also important. Up to 400 filmsper day are routinely screened in busy radiology clinics. The automatedanalysis does not have to be applied on-line; however, it may bedifficult to process large numbers of mammograms overnight if algorithmsare not fast enough. Because the segmentation algorithm has to segmentall candidate structures that may potentially be microcalcifications,its speed is especially relevant. Each film may have several thousandcandidate structures that must be segmented.

The multi-tolerance segmentation algorithm of Shen et al. (L. Shen, etal. “Detection and Classifications of Mammographic Calcifications,”International Journal of Pattern Recognition and ArtificialIntelligence, vol. 7, pp. 1403-1416, 1993), does not use statisticalmodels for local statistics, and its threshold is set automatically.This multi-tolerance, region growing approach uses a growth toleranceparameter that changes in a small range with a step size that depends onthe seed pixel. The structure of interest is segmented multiple timeswith varying tolerance parameters, and in each segmentation, a set ofthree features is computed. The normalized vector differences in thefeature set between successive segmentations are calculated and thesegmentation with minimal difference is selected as the final one.

The active contours model of Kass et al. (Kass, M. et al. “Snakes:Active Contour Models,” International Journal on Computer Vision, pp.321-331, 1988), also provides segmentation without parametricstatistical data models or windows for local statistics, but does relyon several user selected parameters that place some burden on the user.It has been used successfully to determine the boundaries of tissuestructures in data such as ultrasound and MRI images of the heart, andMRI images of the brain, but it has not been applied to the segmentationof microcalcifications. The active contours model starts with an initialcontour placed near the expected boundary and moves the contouriteratively toward the boundary by minimizing an energy function. Thecontour is modeled as a physical flexible object with elasticity andrigidity properties. Its dynamics, dictated by the balance between theseinternal properties and external forces that depend on the image data,satisfy the Euler equations and minimize the corresponding energyfunction. An active contour that is initiated as a closed curve remainsso during iterations and its smoothness can be adjusted by the choice ofparameters.

What is needed is a segmentation method and apparatus withoutstatistical models, local statistics, or thresholds to be selectedmanually, and with significantly lower computational complexity comparedto the multi-tolerance and active contours methods, for enhanced speed.

In particular, what is needed is a method and apparatus to segmentpixels in an image, such as a mammogram, containing a plurality of extradark or extra bright objects just a few pixels in extent, that givesedges similar to those selected by an expert, but does so with fewercomputations and with fewer manually adjustable parameters thanconventional segmentation methods and equipment.

SUMMARY OF THE INVENTION

Therefore it is an object of the present invention to providesegmentation for small features in multidimensional data which definessmall feature edges that correspond closely to those selected by ananalyst but does so with less complexity than the above known methods.

It is another object of the present invention to provide a dataprocessing apparatus that more rapidly provides small feature edges thatcorrespond closely to those selected by an analyst.

It is another object of the present invention to provide computerprogram products that more rapidly provide small feature edges thatcorrespond closely to those selected by an analyst.

It is another object of the invention to identify microcalcifications ina mammogram.

These and other objects and advantages of the present invention areprovided by a method for segmenting a small feature in amultidimensional digital array of intensity values in a data processor.Each small feature includes a local intensity extreme, such as anintensity maximum. An edge metric is computed along each ray of aplurality of multidimensional rays originating at the local intensityextreme. A multidimensional edge point is identified corresponding to amaximum edge metric on each ray. Every point on each ray from the localextreme to the ray edge point is labeled as part of the small feature.The labeling is then spread to an unlabeled point following a hillclimbing procedure requiring that the unlabeled point be adjacent to alabeled point, have a similar or more extreme intensity than the labeledpoint, and be closer than the labeled point to the local extreme.

In another embodiment, the multidimensional array is a digital image,and each point is a pixel. In another embodiment, the digital image is adigitized mammogram and the small feature is a microcalcificationcandidate. In the latter embodiment, microcalcification candidates aresatisfactory segmented in fewer operations than with conventionalsegmentation methods.

In another aspect of the invention, a data processing apparatus segmentsa small feature in a multidimensional digital array of intensity values.The apparatus includes an input for inputting a plurality of intensityvalues arranged along regular increments in each of a plurality ofdimensions and a memory medium for storing the plurality of intensityvalues as a multidimensional digital array. The apparatus includes aprocessor configured to detect a local intensity extreme in themultidimensional digital array, to identify points along a plurality ofrays originating at the local intensity extreme, and to identify one rayedge point on each ray. The ray edge point is associated with a maximumedge metric along the ray. The processor is also configured to label thepoints in the array that are part of the small features. Each point oneach ray from the local intensity extreme to the edge point is labeled,as is an unlabeled point adjacent to a labeled point if the unlabeledpoint has a more extreme intensity than the labeled point and theunlabeled point is closer than the labeled point to the local extreme.Labeling continues until no more unlabeled points can be labeled. Theapparatus also includes an output for providing the labeled points forsubsequent processing.

In another aspect of the invention, a computer program product isprovided for segmenting a small feature in a multidimensional array ofintensities using a computer. The computer program product includescomputer controlling instructions for configuring a computer to computean edge metric along each ray of a plurality of multidimensional raysoriginating at a local intensity extreme. The instructions also identifya ray edge multidimensional point corresponding to a maximum edge metricon each ray. The program also labels every point on each ray from thelocal extreme to the ray edge point, and then labels an unlabeled pointif the unlabeled point is adjacent to a labeled point and the unlabeledpoint has a more extreme intensity than the labeled point, and theunlabeled point is closer than the labeled point to the local extreme.In one embodiment, the instructions are stored in a computer readablememory device. In another embodiment, the instructions are transmittedas electronic signals on a communications line.

The foregoing and other features, aspects and advantages of the presentinvention will become more apparent from the following detaileddescription of the present invention when taken in conjunction with theaccompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

The preferred and example embodiments of the present invention aredescribed with reference to the Drawings in which:

FIG. 1A is a perspective view of the external features of a computerapparatus suitable for one embodiment of the present invention.

FIG. 1B is a block diagram of a computer apparatus that can beconfigured according to one embodiment of the present invention.

FIG. 1C is a perspective view of a sample memory medium for storinginstructions to configure a computer according to another embodiment ofthe present invention.

FIG. 1D is a block diagram of a network that can transmit electronicsignals that configure a computer according to still another embodimentof the present invention.

FIG. 2A is a flow diagram for a method according to yet anotherembodiment of the present invention.

FIG. 2B is a flow diagram following step 270 of FIG. 2A according to afurther embodiment of the present invention.

FIG. 2C is a flow diagram for details of step 260 of FIG. 2A accordingto still another embodiment of the present invention.

FIG. 2D is a flow diagram for an alternative detail for step 260 of FIG.2A according to yet another embodiment of the present invention.

FIG. 3 is a schematic diagram of a local maximum, rays and edges thatresults from steps 210 through 250 of FIG. 2.

FIG. 4 is a schematic diagram of a local maximum, a labeled pixel,adjacent pixels, and a reference line according to one criteria for oneembodiment of step 260 of FIG. 2.

FIG. 5 is a schematic diagram of a local maximum, a labeled pixel, andan adjacent pixel according to a criteria for another embodiment of step260 of FIG. 2.

FIGS. 6A-6D are gray scale photographs showing an actual intensitymaximum as originally provided and then superposed with labeled pixelsafter three stages of the method of FIG. 2 according to the presentinvention.

FIGS. 7A-7D are gray scale photographs showing three actual intensitymaxima as originally provided and then superposed with labeled edgepixels after segmentation based on two conventional methods and thepreferred embodiment of the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

The principles of the present invention will be described next, detailedin terms of preferred and example embodiments with reference to theaccompanying drawings. Whenever possible, the same reference numberswill be used throughout the drawings to refer to the same or like parts.

The explanations of the detailed embodiments are by way of example onlyand are not meant to limit the scope of the invention. The inventionapplies to identifying small structures in any multidimensional array ofregularly spaced intensity values. Here intensity is used in a genericsense representative of measured data values in general, and is notconfined to density of optical energy. Examples of such multidimensionalarrays include gray-scale digital images in which intensity values areregularly spaced in two dimensions, often called rows and columns or yand x, such as the mammogram described in the preferred embodiment. Inthis kind of arrangement, each digital image element is a pictureelement called a pixel. Elevation maps are two dimensional arrays ofheight data, where height is the “intensity.” Other examples ofmultidimensional arrays include color images which can be represented asthree-dimensional arrays of intensity where the third dimension iscolor. Typically, the array would have intensity at only three points inthe color dimension, for example, a red intensity, a blue intensity anda green intensity. Gray-scale video clips can also be consideredthree-dimensional arrays, where each video image frame istwo-dimensional and the third dimension is time. By the same token,color video clips can be considered four-dimensional where the fourdimensions are row, column, color and time. Other examples includemedical imagery where two-dimensional cross sections of a human body areassembled at several positions from head to toe. In this case the thirddimension is height through the subject. By extension, suchthree-dimensional looks can be repeated at uniform intervals of time,making time the fourth dimension. Thus the descriptions that followapply not only to gray scale images of the preferred embodiment, but tomultidimensional arrays of digital data.

A multidimensional point in a multidimensional digital array is locatedby the index of the point in each of the dimensions. Letting D representthe number of dimensions, the location of a multidimensional point P ina multidimensional array can be specified uniquely by a set containing Dindexes as coordinates, {I₁, I₂, I₃, . . . I_(D)}. Where there are onlytwo dimensions, it is common to refer to I₁ as the x coordinate and torefer to I₂ as the y coordinate. There is an implied limit to the numberof allowed positions in each dimension of a finite array. Letting L_(i)represent the maximum number of locations in the i-th dimension of thedigital data array, each index can vary from one to L_(i), inclusive.That is:1≦I_(i)≦L_(i).   (1)The distance, d, between any two multidimensional points, P_(a) andP_(b), with different indices {a₁, a₂, a₃, . . . a_(D)} and {b₁, b₂, b₃,. . . b_(D)}, can be defined as the square root of the sum of thesquares of the differences in their indices. That is,

$\begin{matrix}\begin{matrix}{{d\left( {P_{a},P_{b}} \right)} = {d\left( {{P\left( {a_{1},a_{2},\ldots\mspace{14mu},a_{D}} \right)},{P\left( {b_{1},b_{2},\ldots\mspace{14mu},b_{D}} \right)}} \right)}} \\{= \left. \sqrt{}\left\lbrack {\left( {b_{1} - a_{1}} \right)^{2} + \left( {b_{2} - a_{2}} \right)^{2} + \ldots + \left( {b_{D} - a_{D}} \right)^{2}} \right\rbrack \right.}\end{matrix} & (2)\end{matrix}$

The intensity, f, varies with position in the multidimensional array andmay be represented by the symbol f(P). The intensity f at eachmultidimensional point can be a single value, also called a scalarquantity. Alternatively, the intensity can be a vector of severalvalues, e.g., f(P)={f1(P), f2(P), f3(P)}. For example, the three-colorimage can be treated as a three-dimensional array or can be treated as atwo dimensional image with a three element vector intensity. In thisterminology, the vector elements of the intensity are not used in thecalculation of distance using Equation 2. Instead, the magnitude ofintensity at point P could be any vector magnitude convention such asthe square root of the sum of the squares of the vector components orthe sum of the absolute values of the vector components. Similarly, thedifference in intensity between two points P_(a) and P_(b) would begiven by the magnitude of the difference in the components using anyconventional method.

Thus, though the preferred embodiment is described in which the digitaldata array is an image having two dimensional pixels, each pixel havinga scalar image intensity, the method can readily be extended to multipledimensions using the above relationships. In the following, each pixel Phas a first coordinate represented by x and a second coordinaterepresented by y and an intensity represented by f(P) or f(x,y).Separate pixels are designated by separate subscripts.

Though the invention applies to any imagery, the preferred embodimentssegment two-dimensional images with a gray-scale intensityrepresentative of a mamnmogram. Other two dimensional imagery which thepresent invention can segment include imagery of military scenes inwhich the intensity is responsive to the presence of targets of a firingsystem, such as vehicles to be fired upon by a missile.

The invention is related to finding small objects in a multidimensionalarray. In this context small means objects affecting the intensity inseveral points in one dimension of the array but not many thousands ofpoints in each dimension. Other, statistical and textural segmentationprocedures, for example, are expected to be more useful as the number ofpoints in a feature increases. It is characteristic ofmicrocalcifications in mammograms and distant targets in militaryscenarios that only several pixels are contained in the object to besegmented. It is also anticipated that many other features to bedetected in radiographs and sonograms of biological bodies also involveonly several pixels. The present invention is expected to performespecially well for these applications.

The methods and procedures discussed herein are intended to be performedby data processing systems or other machines. Though described in termsthat can be interpreted to be performed by a human operator, suchperformance is neither required nor likely to be desirable. Multipletedious computations with high accuracy are required that are unsuitablefor practical implementation by human beings. Also, the invention can beimplemented in computer or other hardware, the structure of which isevident from the following descriptions.

Also herein, the procedures will be described as the manipulation ofvalues, symbols, characters, numbers, or other such terms. Though suchterms can refer to mental abstractions, herein they are used asconvenient expressions for physical signals such as controllablechemical, biological, and electronic and other physical states that canbe used to represent the values, symbols, characters, numbers, or othersuch terms.

FIG. 1A illustrates a computer of a type suitable for carrying out theinvention. Viewed externally in FIG. 1A, a computer system has a centralprocessing unit 100 having disk drives 110A and 110B. Disk driveindications 110A and 110B are merely symbolic of a number of disk drivesthat might be accommodated by the computer system. Typically these wouldinclude a floppy disk drive such as 110A, a hard disk drive (not shownexternally) and a CD-ROM drive indicated by slot 110B. The number andtype of drives vary, typically, with different computer configurations.The computer has a display 120 upon which information is displayed. Akeyboard 130 and mouse 140 are typically also available as inputdevices.

FIG. 1B illustrates a block diagram of the internal hardware of thecomputer of FIG. 1A. A bus 150 serves as the main information highwayinterconnecting the other components to the computer. CPU 155 is thecentral processing unit of the system, performing calculations and logicoperations required to execute programs. Read-Only-Memory 160 andRandom-Access-Memory 165 constitute the main memory of the computer.Disk controller 170 interfaces one or more disk drives to the system bus150. These disk drives may be floppy disks drives, such as 173, internalor external hard drives, such as 172, or CD-ROM or DVD (digital videodisk) drives such as 171. A display interface 125 interfaces a display120 and permits information from the bus to be viewed on the display120. Communications with external devices can occur over communicationsport 175.

FIG. 1C illustrates an exemplary memory medium which can be used withdrives such as 173 in FIG. 1B or 110A in FIG. 1A. Typically, memorymedia such as a floppy disk, or CD-ROM, or DVD, will contain the programinformation for controlling the computer to enable the computer toperform its functions in accordance with the invention.

FIG. 1D is a block diagram of a network architecture suitable forcarrying data and programs over communication lines in accordance withsome aspects of the inventions. A network 190 serves to connect a usercomputer or client computer 110 with one or more servers such as server195 for the download of program and data information. A second user on asecond client computer 100′ can also connect to the network via anetwork service provider, such as ISP 180.

In general, small objects in images may have an intensity level that iseither lower or higher than a surrounding background. An intensitymaximum with levels higher than the background is called a localmaximum, and an intensity minimum with intensity levels below thebackground is called a local minimum. Both maximum and minimum areencompassed by the term intensity extreme. Thus, in general, the targetobjects in an image or multi-dimensional array encompass intensityextremes. Both are capable of being segmented according to the presentinvention. For the sake of serving as an example, the followingdescription generally considers the preferred embodiment in whichmicrocalcifications are evident as local maxima in intensity, and themethod will be called a hill climbing method; however, segmenting alocal minimum is also anticipated using the hill climbing method. In thefollowing discussion, when a first point has an intensity equaling theintensity of the local extreme or between the intensity of the localextreme and the intensity of a second point, the first point is said tohave a more extreme intensity than the second point.

FIG. 2A shows the method according to one embodiment of the presentinvention. A local brightness maximum, characteristic of amicrocalcification, is identified at pixel P₀ in an image at step 210.Next, a plurality of rays is defined that emanate from that localmaximum pixel P₀ as illustrated in step 220. FIG. 3 illustrates fivesample rays 320 emanating from a local maximum 310. Referring again toFIG. 2A, an edge metric is computed for each pixel along each ray instep 230. Then in step 240, a ray edge pixel on the ray is identifiedbased on a maximum edge metric. Then the pixels on the ray from thelocal maximum to the ray edge pixel, inclusive, are labeled as belongingto the object or feature in step 250. Additional pixels belonging to thefeature are labeled if they are adjacent to a labeled pixel and if theunlabeled pixel satisfies intensity and distance criteria describedlater. These criteria implement the unique hill climbing procedure ofthe present invention. This growth of labeled pixels is indicated bystep 216 260. In step 270, every unlabeled pixel next to a labeled pointis examined using the criterion in step 260 until no further points canbe labeled.

FIG. 2B shows steps that follow step 270 in another embodiment of thepresent invention. Here each of the labeled pixels is checked in step275 and those labeled pixels adjacent to an unlabeled pixel arerelabeled as an edge pixel of the small feature. This completes thelabeling associated with one of the small features in the image; and, instep 280, control is returned to step 210 until no local maximum remainsunlabeled or unsegmented in the image. In yet another embodiment of theinvention, small features identified in the image can be joined in step285 if those pixels are within a joint distance. Additional detailregarding the steps shown in FIGS. 2A and 2B are provided with referenceto FIGS. 2C through 5.

According to the present invention, the segmentation is based on theexperience that, in a given array, the edge of a small feature to besegmented is a closed contour around a local intensity extreme pixel P₀.In the preferred embodiment, the local intensity extreme is selected asthe pixel with an extreme intensity (maximum or minimum) in a region thesize of the expected size of the small feature or object. The regionshould have the same number of dimensions as the data array, just fewerpixels. In other words, the region is defined as a sub-array of themultidimensional size equal to the expected size of the feature. In thecase of mammograms, this sub-array is a square that is about 100 pixelsin x and 100 pixels in y when the resolution of the image is about 25microns per pixel. To avoid selecting local extremes that areinsignificant, the extreme is also required to achieve a certainabsolute value—above a pre-set bright threshold in the case of a localmaximum, or below a pre-set dark threshold in the case of a localminimum.

A pixel P on a ray is considered to be on the edge of a small object ifit provides a maximum edge metric in a line search on a ray originatingfrom the local extreme pixel and moving in a direction k. The edgemetric may be defined as the change in intensity with each succeedingpixel in the direction k or by a Sobel operator centered on the pixel,or by any known edge metric. However, in the preferred embodiment with alocal maximum, a ray edge pixel is found that more closely correspondsto that selected by expert analysis when the edge metric is a slopedefined according to equation 3.

$\begin{matrix}{{S(P)} = \frac{{f\left( P_{0} \right)} - {f(P)}}{d\left( {P_{0},P} \right)}} & (3)\end{matrix}$For each pixel P around this local maximum P₀ the slope has a value S(P)where f(P₀) is the intensity, e.g., the gray scale value, at the localmaximum pixel P₀, and f(P) is the intensity at pixel P, and d(P₀, P) isthe distance between the local maximum pixel P₀ and the pixel P. Ingeneral, to extend to the case where P₀ is a local minimum, the absolutevalue of the numerator is used. The notation d(P₁,P₂) here indicates theabsolute value of the distance between two points P₁ and P₂. Let P_(n)represent the nth pixel along a ray in a direction k. The n varies from0 at the local maximum to N−1 at the Nth consecutive pixel along theray. The number N is not a critical choice as long as it is larger thanthe number of pixels expected to lie between the local maximum and theedge of the largest structures of interest. Referring to FIG. 3, Nshould be the number of pixels extending half the length of the arrow330 indicating the maximum expected size of a small feature, forexample. Among the pixels P_(n), the pixel at which S(P_(n)) is maximalis considered to be an edge point in that direction and is denoted bye(k). In the preferred embodiment, the ray search is applied in manyequally spaced directions originating from the local maximum pixel,resulting in a set of ray edge pixels e(k) where k varies from 1 to K,the number of directions for which rays are computed. In the preferredembodiment, as shown in FIG. 3, K equals 16. For each direction k, theedge pixel and all pixels between the local maximum and the edge pixele(k) are labeled as belonging to the object associated with the localmaximum pixel P₀. This results in K radial lines of labeled pixels 350,as shown in FIG. 3. These labeled pixels are used as seeds or referencepixels for growing a region to identify all the pixels of the object.

To identify all pixels lying within a contour including the edge pointse(k), the region should grow essentially on pixels with more extremeintensity (e.g., increasing intensity) and toward the local extreme(e.g., local maximum). From any labeled pixel taken as a referencepoint, the region can grow to an adjacent unlabeled pixel if this newpixel satisfies some particular conditions. In the case of data arrayswith more than two dimensions, adjacent points to a labeled point arethose whose indices are all within one of the corresponding indices ofthe labeled point. Referring to FIG. 4, the reference pixel is thelabeled pixel 420 and the eight adjacent pixels are numbered clockwisefrom the diagonally upper left pixel as pixel 1 through 8. These eightpixels are considered eight-connected with the labeled pixel 420. Asubset of these adjacent pixels is the four-connected set of pixels towhich pixels labeled 2, 4, 6 and 8 belong. With respect to the referenceor labeled pixel 420, an eight-connected adjacent or neighbor pixel ischecked. If the neighbor pixel is already labeled, it has already beendetermined that the neighbor pixel is on the object. If the neighborpixel P is not labeled, then it has to satisfy the following conditionsto be labeled.

IF f(P)≧f(P_(r)) then P must be in a position that constitutes a stepfrom P_(r) toward P₀.

IF f(P)<f(P_(r)), then P should be closer to P₀ than P_(r) is to P₀ bymore than a minimum distance called an inclusion tolerance distance.

All pixels labeled during the process are used as reference pixels. Themethod stops when no pixel can be appended as shown in step 270 of FIG.2A. The step for labeling unlabeled pixels is illustrated in FIG. 2A asstep 260.

The intensity and distance criterion referred to in step 260 are nowdescribed with reference to FIGS. 2C and 2D, which each show one of thetwo alternative criteria used in the present hill climbing method andapparatus. In both these figures, the first condition checked is theintensity f(P) of the unlabeled point P compared to the intensityf(P_(r)) at the reference pixel P_(r), as shown in step 262.

Most microcalcifications have an intensity that decreases monotonicallyfrom the local maximum toward the edges. However, in some cases, thismay not be true, and the growth toward the local maximum may need toinclude new pixels that have lower values or less extreme values thantheir labeled referenced pixels. As long as this is done strictly towardthe local extreme, growth in an unwanted direction is avoided. That is,if the unlabeled pixel P is much closer to the local maximum (orminimum) than is the labeled referenced pixel P_(r), then the unlabeledpixel P is considered engulfed by the object and is labeled even if itsintensity f(P) is less extreme than f(P_(r)). The distance by which theunlabeled point must be closer than the labeled point to be engulfed bythe object is called the inclusion tolerance distance. In this and thefollowing discussions, the difference in distances between the labeledand unlabeled points to the local maximum P₀ is represented by G givenin Equation 4.G=d(P₀, P)−d(P₀, P_(r))   (4)When the unlabeled pixel P is closer to the local maximum P₀ than theunlabeled pixel P_(r), then G is negative. Therefore, the negative of Gis compared to the inclusion tolerance to determine if the unlabeledpixel is close enough to the local extreme to be engulfed, as shown instep 263 of FIGS. 2C and 2D. In the preferred embodiment, the inclusiontolerance is one pixel. Thus, lower intensity pixels closer to the localmaximum than the already labeled point P_(r) by more than one pixel areclose enough to be labeled. That is, a new pixel P with intensity f(P)less extreme than the intensity f(P_(r)) of the referenced pixel P_(r)is appended to the region if its distance to the local extreme is suchthat −G is ≧ the inclusion tolerance distance, as shown in step 265 ofFIGS. 2C and 2D. If the unlabeled pixel with less extreme value is lessthan the inclusion tolerance closer to the local extreme or is fartherfrom the local extreme, then the unlabeled pixel is not labeled, asshown in step 265 of FIGS. 2C and 2D.

The other branch from step 262 in FIGS. 2C and 2D is followed when theadjacent pixel P that is unlabeled has an intensity that is greater thanor equal to the intensity of the labeled pixel P_(r). This correspondsto the condition in the case of a local minimum that the unlabeled pixelhas a lower intensity than the labeled pixel P_(r). That is, the “yes”branch is followed from box 267, in general, if the unlabeled pixel Phas an intensity that is no less extreme than the intensity at thelabeled pixel P_(r). Each of two different criteria can be used todetermine whether the unlabeled pixel P is in a position thatconstitutes a step from the labeled pixel P_(r) toward the extreme pixelP₀.

The first criterion, Criterion 1, is indicated in FIG. 2C and step 264aand is based on the angle of the line perpendicular to the line segmentconnecting the local extreme P₀ with the reference pixel P_(r). The lineperpendicular to the segment connecting the local extreme to the labeledpixel is called the reference line 430 and is shown in FIG. 4. Forarrays of more than two dimensions, the reference would be a surfacewith a number of dimensions at least one dimension less than themultidimensional array. The numbered pixels of FIG. 4 are approved forappending to the small feature if they fall within the list of approvedpixels listed in Table 1 for the quadrant in which the angle θ variesfrom 0-90°. The angle θ

TABLE 1 Criterion 1 for First Quadrant. xr yr 0 Approved Pixels x_(r) =x_(o) y_(r) < y_(o) 0° 1, 2, 3, 4, 8 x_(r) > x_(o) ″ 0 < tan θ ≦ ⅓ 1, 2,3, 4, 8 ″ ″ ½ < tan θ < 1 1, 2, 3, 8 ″ ″ tan θ = 1 1, 2, 3, 7, 8 ″ ″ 1 <tan θ ≦ 3 1, 2, 7, 8 ″ ″ 3 < tan θ < ∞ 1, 2, 6, 7, 8 ″ y_(r) = y_(o) 90°1, 2, 6, 7, 8between the reference line 430 and the x-axis is also shown in FIG. 4.The first two columns of Table 1 show the relationship between thecoordinates of the reference pixel x_(r) and y_(r) of P_(r) and theirrelationship to the coordinates x₀ and y₀ of the local maximum P₀. Fordifferent values of the angle θ or its tangent, tan θ, different of thenumbered pixels in FIG. 4 are approved. Table 1 captures the conditionthat the unlabeled pixel P and the local maximum P₀ must lie on the sameside of the reference line 430. Among the eight pixels that surround areference pixel, only some will meet the spatial criterion of Criterion1, depending on the angel θ of the reference line. The angle θ ismeasured positive counterclockwise from the x-axis. The allowable pixelsfor values of θ in the other three quadrants are obtained in asymmetrical manner. An extended table would have to be drafted for dataarrays of greater than two dimensions.

AsReferring to FIG. 5, as an alternative for the constraint C1 describedabove and summarized in Table 1, Constraint 2 can be used to determinewhether a neighboring pixel should be labeled. Constraint 2 is morereadily extensible to more than two dimensions. Referring to Equation 4defining the distance difference G, most allowable pixels described byCriterion 1 yield a negative G value. However, some pixels generate apositive G value. These positive G pixels are the pixels that provide astep, from the reference pixel P_(r), approximately parallel to thereference line. This type of growth through pixels is especiallydesirable around the edge of the small structure. The largest values ofG are associated with diagonal pixels and occur at the edge of thesmallest features to be segmented. Furthermore, among all possible pixelconfigurations, the value of G is maximal when the reference line angleθ is 45° or 135° and the new pixel P is diagonally connected to thereferenced pixel P_(r). This maximal value is also obtained for otherhomologous arrangements of the three pixels. A positive threshold G_(t)for G can be used instead of Criterion 1. Consider an approximatelycircular object 2N pixels wide. On the edge of such an object, thehighest value for G, called G_(max), will equal (√(N²+2))−N. The smallerN, the larger G_(max) will be. An appropriate threshold for G can be setby using the width of the smallest object of interest. Therefore, analternative way of constraining the expansion of pixels away from thelocal extreme is to allow only new pixels which provide a value of G ofat most G_(max). That is, set G_(t)=G_(max). This threshold, G_(t), canbe considered an expansive tolerance distance. Criterion 2 can be statedas: G must be less than or equal to the expansive tolerance distance,G_(t). For example, mammograms with pixels of 25 microns andmicrocalcification candidates having structures as small as 0.25 mmacross, yield N=5; so, G_(t)=G_(max)=0.196.

The preferred embodiment determines 16 ray edge pixels around theobject, and segments with the hill climbing procedure described. Asindicated in step 270 of FIG. 2A, each appended pixel is labeled and isused as a reference pixel itself during growth. The growth stops when nopixel can be appended. Once no more new pixels can be labeled, eachlabeled pixel is examined to identify edge pixels of the small featurein step 275 in FIG. 2B. The edge pixels of the small feature aredetermined to be all labeled pixels that are four-connected to anunlabeled pixel after no further pixels can be added.

After every object has been segmented and its outer edge pixels defined,larger features may be discernable. The larger features can beconstructed where the small features abut or overlap slightly. The stepof joining small features together into a larger feature is depicted instep 285 of FIG. 2B. Depending on the larger feature being assembled,the criterion for joining small features can be that the small featuresshare edge pixels, or that the edges overlap so that the edge of onesmall feature is an interior labeled pixel of another small feature. Itis also possible that features be joined that do not touch or overlap,provided they are sufficiently close together. A tolerance called a joindistance can be used to determine how close the edges should be to eachother in order to combine the small features into one or more largerfeatures. In this case, all small features are joined where the edgepixels of two different small features are within a join distance.Overlapping pixels are covered by this criterion as are features whoseedge pixels coincide. By setting the joined distance to 0 the edgecoincidence is required; and by setting the join distance negative,overlapping can be required.

EXAMPLES

To determine whether the results of the present invention provide edgesof small features that are useful in interpreting mammograms and indoing so with fewer computations than other methods, several experimentswere performed with actual mammograms. The correctness of the edgedetermined by the present invention is measured by its similarity to theedges determined by an analyst, and its ability to discriminate amongthe candidate microcalcifications in subsequent processing. Otheradvantages of the preferred embodiment are measured using the complexityor number of computations involved in the procedure, and the timerequired to execute the procedure on a computer.

Example 1

Five mammograms containing subtle microcalcification clusters were usedto evaluate the algorithms for data that would warrant the use of anautomated system. Mammograms without magnification were used; and thebreast images covered an area that ranged between 12 cm×6 cm and 21cm×11 cm. The location of individual microcalcifications was indicatedby an experienced mammographer. These 5 mammograms contained 15 clusterswith a total of 124 microcalcifications, yielding about 8microcalcification per cluster. The number of microcalcifications percluster ranged between 3 and 18. The size of microcalcifications rangedbetween 0.25 mm and 1 mm wide, with more than 90% being smaller than 0.5mm. Mammograms were digitized with a Howtek D4000 drum scanner using aspatial resolution of 25 microns per pixel and 12-bit A/D conversion,with an optical dynamic range of 0-3.5 optical depths (O.D.).

The multi-tolerance region growing procedure grows a region around aseed pixel by appending 4-connected pixels P that satisfy:(1+τ)(F_(max)+F_(min))/2≧P≧(1−τ)(F_(max)+F_(min))/2   (5)where τ is the tolerance parameter, and F_(max) and F_(min) are thecurrent maximum and minimum values in the region grown that far. Thevalue of τ is not manually selected by the user; the best value isautomatically determined for each segmented structure by repeating thegrowth with multiple values of τ between 0.01 and 0.4 with steps ofs=1/v, where v is the 8-bit value of the seed pixel. Three features areextracted from each region grown with a different tolerance level: shapecompactness, center of gravity, and size. The algorithm determines thevalue of T that results in the minimal change in the vector of thesethree features with respect to the previous τ value in the sequence bycomputing a normalized difference between consecutive vectors. Thevector with minimal difference indicates the best choice of τ.

The segmentation outcome of the multi-tolerance region growing procedureon 5 subtle microcalcification candidates depended partly on theintensity structure of the microcalcification. When the intensitytransition from the edge to the background was relatively abrupt, thesegmented region coincided closely to the visually perceived edge. Whenthe intensity at the edge decreased gradually toward the backgroundlevel, this algorithm generally produced a relatively large region.Nevertheless, the growth was consistently contained, i.e. it did notgrow to an unacceptable size and it generated boundaries that can beused as an estimate of the immediate background around themicrocalcification.

The active contours model represents the contour points asv(s)=(x(s),y(s)) The contour is obtained by minimizing the energyfunctional:E[v(s)]=∫_(Ω)Eint[v(s)]+PE[v(s)]+Eext[v(s)]ds   (6)where E_(int) is the internal energy due to the elasticity and therigidity, PE is the potential energy obtained from the image data,E_(ext) is the energy of external forces that can be applied to thecontour. The integration is performed over the entire contour Ω. Theinternal energy is expressed by:E_(int)=w₁|v′(s)|²+w₂|v″(s)|²   (7)where w₁ and w₂ are coefficients that control the elasticity andrigidity, respectively, and primes denote differentiation. The choice ofpotential energy depends on the application; it is typically thenegative squared gradient magnitude, and is so used for mammograms.

The active contour that minimizes E(v) satisfies the Euler-Lagrangeequation:−(w₁v′)′+(w₂v″)″=F(v)   (8)where F(v) represents the force due to the combined effects of thepotential energy and external energy. In this study we implemented theballoon forces and the image force normalization suggested, resulting in

$\begin{matrix}{{F(v)} = {{k_{1}{n(s)}} - {k_{2}\frac{\nabla{PE}}{{\nabla{PE}}}}}} & (9)\end{matrix}$where n(s) is the unit vector normal to the contour at point v(s),oriented toward the outside of the contour, k₁ is the magnitude of theballoon inflation force, and k₂ is the coefficient of the normalizedimage force. The value of k₂ is selected to be slightly larger than k₁to allow edge points to stop the inflation force.

The numerical solution was implemented using finite differences and theiterative evolution as suggested:(I+τA)v_(t)=(v_(t−1)+τF(v_(t−1)))   (10)where I is the identity matrix, τ is the time step, A is thepentadiagonal matrix obtained with the finite difference formulation ofE_(int), v_(t) is the active contour vector at time t, and F(v_(t)) isthe external force vector at time t. We used the negative squaredmagnitude of the image gradient as the potential energy. Pixels detectedwith an edge detector were not used in this study. The gradient of theimage was computed with the Sobel operator.

The initial position of the contour was set automatically for eachstructure to be segmented. Since each structure of interest is a localintensity extreme, pixels were selected that were local maxima acrossthe entire image. Each local maximum was used to segment a region aroundit. The width of the smallest microcalcifications considered in thisstudy was about 0.25 mm and the majority of the microcalcifications inour database had widths in the range 0.3 to 0.5 mm. A circle of 0.2 mmdiameter around the local maximum pixel was used as the initial positionof the active contour. The initial contour points were 248-connectedpixels forming this circle.

The selection of parameters for the active contour segmentation requiredsome trial and error to obtain good segmentation. The segmentation ofthe same 5 subtle microcalcification candidates was performed usingdifferent active contours parameters. First, following therecommendations of Cohen (Cohen, L. D. “On Active Contour Models andBalloons,” CV GIP: Image Understanding, vol. 53, pp. 211-218, 1991), weselected the values of w₁ and w₂ as a function of the spatialdiscretization step size h, such that w₁ was of the order of h² and w₂was of the order of h⁴(w₁=6, w₂=40). Then τ was also set to 0.1. When k₁and k₂ were relatively small (2 and 4), the image force and the balloonforce did not act sufficiently on the active contour, producing contoursthat were only slightly different than the initial position. When thesetwo parameters were increased (14 and 16), the resulting segmentationwas very close to that expected visually. Increasing these parametersfurther (24 and 26) increased the combined effect of image gradient andballoon forces, producing contours that extended beyond the expectededges. Within this range, segmentation with the active contour model wasnot very sensitive to the values of the other parameters. The effect ofdoubling w₁ to 12, is that contours became slightly smaller due to theincreased stiffness of the active contour model. Sensitivity to w₂ wasalso low. When w₂ was doubled to 80, the contours became slightlysmoother due to the increased rigidity of the model.

The segmentation steps of the hill climbing approach of the presentinvention are illustrated in FIG. 6. FIG. 6A shows a microcalcificationcandidate that has a width of about 0.3 mm. The 16 ray edge points 624determined by the radial line search of the hill climbing algorithm areshown in FIG. 6B. The region grown using spatial Constraint 1 is in FIG.6C. The region grown with spatial Constraint 2 was identical for thismicrocalcification candidate. The edge pixels 642 of the entiremicrocalcification candidate are shown in FIG. 6D. The segmentation ofmicrocalcifications by the hill climbing method produces outcomes usingthe spatial Constraints 1 and 2 that were almost identical. In thisstudy, about a quarter of microcalcifications were segmented identicallyby the two spatial constraints and the rest differed by a few pixels,resulting in a negligible change over the entire microcalcification.Both spatial constraints directed the growth of the regionssuccessfully, resulting in regions that were compatible with visualinterpretation.

The differences between the three methods are illustrated in FIG. 7.Three subtle microcalcifications candidates are shown in FIG. 7A. Whenthe contrast of a microcalcification candidate was relatively low, orparts of it exhibited a very gradual decrease in intensity toward thebackground, the multi-tolerance algorithm (FIG. 7B) segmented a largerregion than those of the other two algorithms. Good segmentation withactive contours (FIG. 7C) was obtained using w₁=6, w₂=40, τ=0.1, k₁=14and k₂=16, for all microcalcifications candidates of this study. Usingthese parameters, segmentation with active contours provided edges 735that were smoother than edges 725 and 745 produced by segmentation withthe other two methods. The selection of w₁ and w₂ provided theflexibility needed to adapt relatively well to the shape of diversemicrocalcifications candidates. The elasticity level allowed the contourto grow to the highest gradient locations when the segmented structureswere relatively large, and the rigidity level allowed the contour todevelop sharp bends dictated by the data in some microcalcifications.The edges 745 of regions grown by the hill climbing algorithm shown inFIG. 7D were not as smooth as those 735 of the active contours, but theconvolutions were consistent with visually perceived edges aroundmicrocalcifications candidates.

Example 2

Segmentation of microcalcification candidates serves as an initial stepfor discriminating between the population of microcalcifications andthat of background structures. The discrimination potential of eachsegmentation algorithm was quantified using features extracted fromstructures segmented around all the local maxima in the 5 mammograms.These structures consisted of the 124 microcalcifications mentionedabove and 2,212 background structures segmented in the same mammograms.Four characteristics were used to assess the discrimination potential inthis study.

1. Contrast was measured as the gray level difference between the localmaximum pixel P₀ in the structure, and the mean of pixels around itsedge.

2. Relative contrast was computed as the ratio of the contrast to thevalue at the local maximum.

3. Area was computed as the number of labeled pixels in the grownregion.

4. Edge sharpness was the mean of the gradient computed with a Sobeloperator across all edge pixels. The Sobel operator is a mask whichweights the eight neighbors of a pixel to compute a sum proportional tothe gradient x, or the y gradient, or total gradient.

The discrimination ability of each feature was determined separatelyusing the area under a receiver operating characteristic (ROC) curveobtained with that feature. The ROC curve pots the percentage ofcorrectly detected microcalcifications against the percentage ofdetected background structures as a detection threshold is changed. TheROC curve area is higher when the feature has distributions that aremore separable for a given property. When both populations overlapcompletely, the ROC curve area is 0.5. In general, effectivediscrimination power is indicated by a value above 0.8. Table 2summarizes the results for all three procedures. The area feature hadvery low discrimination power for all three algorithms, indicating thatthe two types of structures cannot be discriminated well on the basis oftheir area segmented. However, the other

TABLE 2 Multi-tolerance Active Hill Region Growing Contours ClimbingContrast 0.80 0.82 0.83 Relative Contrast 0.83 0.90 0.90 Area 0.63 0.600.54 Sharpness 0.80 0.85 0.85three features suggested good discrimination potential for all threealgorithms. A comparison among algorithms shows that both the hillclimbing method of the present invention and the active contoursalgorithm provide segmentation with the same discrimination power, andthey both perform slightly better than the multi-tolerance segmentation.Thus, the hill climbing method produces edges as good as the bestproduced by the conventional approaches tested.

The significant advantage of the hill climbing algorithm is its speed.While the multi-tolerance algorithm provides a good solution to avoidthe use of statistical models, local statistics estimators and themanual selection of threshold, its cost is multiple segmentations of thesame structure and computation of features during the segmentation ofeach structure. Furthermore, in some cases, this algorithm segmentsregions that are somewhat larger than expected. Consequently, the timerequired for segmentation of a mammogram with this algorithm is high.The segmented regions were comparable to those of the other twoalgorithms in many cases. The differences were caused by the fact thatthe growth mechanism of this algorithm is constrained only by anintensity range criterion applied to a new pixel. In contrast, activecontours are constrained by internal forces that regulate the growthaway from the local maximum, and hill climbing has an inward growthmechanism based on edge points.

The active contours also circumvent the statistical and manual thresholdselection issues for each mammogram, but the selection of theoperational parameters for a set of mammograms requires some trial anderror. However, when an appropriate set of parameters is determined, itappears to be valid for a wide range of microcalcifications so it neednot be modified with each mammogram. The choice of negative squaredgradient magnitude as the image energy function seems to be anappropriate one to segment microcalcifications.

Example 3

The computational complexity c_(m) of the multi-tolerance region growingalgorithm is of the order O(4smo) where s is the number of steps in thetolerance search, m is the number of pixels in the region, and o is thenumber of operations per pixel. The factor 4 is included because thealgorithm visits the 4-connected neighbors for each pixel in the region.Considering 125 to be an average intensity value for the local maximum,the average step size is 0.008 resulting on the average in about s=50steps to cover the range 0.01 to 0.4. The average size of segmentedstructures is about 200 pixels. At each pixel the computations performedinclude intensity comparisons, update of F_(max) and F_(min), andcalculation of the center of gravity. Considering about 12 operationsper pixel on the average, the numerical estimate for the average numberof operations per segmentations is c_(m)=480,000.

The computational complexity c_(a) of the active contour model isO[2(n+n²)t] where n is the number of contour points, and t is the numberof iterations. The factor of 2 is included due to the fact that the xand y coordinates of each contour point are computed separately, withidentical operations. At each iteration, order n computations are neededto determine the normal vectors, and order 2n² operations are needed toperform a matrix multiplication. In this study 24 contour points wereused, and the number of iterations depended on the size of thestructure. On the average however, the active contour model converged inabout 20 iterations. This resulted in an average value of c_(a)=47,040,a factor of ten improvement over the multi-tolerance method.

The complexity c_(h) of the hill climbing method is O(KN+8 m) where K isthe number of radial directions from the local maximum, N is the numberof pixels searched in each direction, and m is the number of pixels inthe grown region. A factor of 8 is included since all 8 neighbors ofeach pixel are visited. In this study K was 16 and N was 40, andconsidering an average structure size of m=200, the average estimate ofthe number of operations is c_(h)=2,240, a factor of 20 improvement overthe active contour methods, and 200 over the multi-tolerance method. Theproportions of c_(m), c_(a) and c_(h) are approximately 214:21:1respectively, with hill climbing far less complex than the other twomethods.

Example 4

The speed of the different methods was compared using a section of amammogram containing 456 local maxima, 35 of which were inmicrocalcifications. The sizes of microcalcifications ranged between0.25 mm and 0.5 mm. The times to complete the segmentation of thissection of mammogram using the three algorithms implemented in C on a 10million floating point operations per second (MFLOPS), IBM 6000 computerwere 17 minutes 47 seconds for the multi-tolerance algorithm, 1 minute47 seconds for the active contours, 7 seconds for hill climbing withspatial Constraint 1, and 5.4 seconds for hill climbing with spatialconstraint 2.

Hill climbing with spatial Constraints 1 and 2 yielded practicallyidentical segmentations; but the method was about 20% faster usingspatial constraint 2, resulting in 11.8 ms on average for segmenting astructure, as opposed to 15.3 ms obtained with spatial Constraint 1.

A common technique to determine the edges of an object uses an edgeenhancement algorithm such as the Sobel operator, thresholding toseparate the pixels on edges, and pixel linking to string edge pixelsthat belong to the same object. Selection of the threshold is critical,and linking poses problems in segmenting microcalcifications becausethere are many closely spaced small structures in a background that arelikely to produce considerable numbers of edge pixels. The hill climbingmethod of the preferred embodiment determines edge points that are onthe edge of the same object by virtue of the radial line searchemanating from the same local maximum. It does not require a thresholdto separate edge pixels because the slope in Equation 3 is referred tothe local maximum and is greatest at pixels that are on, or very near,the visually perceived edges. Finally, the hill climbing method avoidssome pitfalls of the region growing mechanism by growing a regioninward, toward the local maximum.

There has been disclosed a segmentation method and apparatus for dataarranged in a multidimensional array which overcomes the problems of theprior art. Although the present invention has been described above byway of detailed embodiments thereof, it is clearly understood thatvariations and modifications may be made by one of ordinary skill in theart and still lie within the spirit and scope of the invention asdefined by the appended claims and their equivalents.

1. A method for segmenting a small feature in a multidimensional digitalarray of intensity values in a data processor, the method comprising:computing an edge metric along each ray of a plurality ofmultidimensional rays originating at a local intensity extreme;identifying a multidimensional edge point corresponding to a maximumedge metric on each said ray; labeling every point on each said ray fromsaid local intensity extreme to said edge point; and labeling anunlabeled point if the unlabeled point is adjacent to a labeled pointand the unlabeled paint point has a more extreme intensity than thelabeled point and the unlabeled point is closer than the labeled pointto the local intensity extreme.
 2. The method of claim 1 whereinintensity is a vector of values and an edge metric is a magnitude of avector difference in intensities between two points along each said raydivided by a multidimensional distance between the same two points. 3.The method of claim 1 further comprising additionally labeling anunlabeled point if the unlabeled point is adjacent to a labeled pointand has a more extreme intensity than the labeled point and is nofarther from the local intensity extreme than the sum of a distance fromthe labeled point to the local intensity extreme plus an expansivetolerance distance less than the spacing between adjacent points.
 4. Themethod of claim 1 further comprising also labeling an unlabeled point ifthe unlabeled point is adjacent to a labeled point and the unlabeledpoint has a less extreme intensity than the labeled point and theunlabeled point is closer than the labeled point to the local intensityextreme by an inclusion tolerance distance.
 5. The method of claim 4,wherein the inclusion tolerance distance is about a spacing distancebetween adjacent points in the array or more.
 6. The method of claim 1,wherein the edge metric at a ray point along each ray is computed as thequotient of the absolute value of an intensity difference between thelocal intensity extreme and the ray point divided by the absolute valueof a distance between the ray point and the local intensity extreme. 7.The method of claim 1, wherein a ray length of each said ray is scaledby an expected size of a small feature.
 8. The method of claim 1,wherein the local intensity extreme is a point with the maximumintensity among a subarray of the multidimensional digital array ofintensity values, the subarray having a certain multidimensional size,and the intensity of the local intensity extreme exceeds a brightthreshold intensity.
 9. The method of claim 8, wherein the certainmultidimensional size is an expected size of a small feature.
 10. Themethod of claim 1, wherein the local intensity extreme is a point withthe minimum intensity among a subarray of the multidimensional digitalarray of intensity values, the subarray having a certainmultidimensional size, and the intensity of the local intensity extremeis less than a dark threshold intensity.
 11. The method of claim 10,wherein the certain multidimensional size is an expected size of a smallfeature.
 12. The method of claim 1, wherein the multidimensional arrayis a digital image, and each point is a pixel.
 13. The method of claim12, wherein the digital image is a digitized mammogram and the smallfeature is a microcalcification candidate.
 14. The method of claim 12,wherein the digital image is a video frame of a military scene and thesmall feature is a candidate target of a tiring firing system.
 15. Themethod of claim 1, wherein said labeling continues until no furtherunlabeled point can be labeled.
 16. The method of claim 15, furthercomprising relabeling a labeled point as a feature edge point if anadjacent point is an unlabeled point.
 17. The method of claim 16,further comprising joining a plurality of small features into acomposite feature when a feature edge point from one small feature ofthe plurality of small features is within a join distance of a featureedge point of another small feature of the plurality of small features.18. A method for segmenting a small feature in a multidimensionaldigital array of intensity values in a dataprocessor, the methodcomprising: computing an edge metric along each ray of plurality ofmultidimensional rays originating at a local intensity extreme:;identifying a multidimensional edge point corresponding to a maximumedge metric on each said ray:; labeling every point on each said rayfrom said local intensity extreme to said edge point; labeling anunlabeled point if the unlabeled point is adjacent to a Labeled labeledpoint and the unlabeled point has a more extreme intensity than thelabeled point and the unlabeled point is closer than the labeled pointto the local intensity extreme:; and additionally labeling an unlabeledpoint if the unlabeled point is adjacent to a labeled point and has amore extreme intensity than the labeled point and is no farther from thelocal intensity extreme than the sum of a distance from the labeledpoint to the local intensity extreme plus an expansive tolerancedistance less than the spacing between adjacent points; wherein anexpected size of a small feature is twice an integral number N times aspacing distance between adjacent points in the array, N is greater than1, the maximum value of the difference in distances between the labeledpoint and the unlabeled point to the local intensity extreme(Gmax)=−N+√{square root over (N²+2))}, and the expansive tolerancedistance is less than about Gmax.
 19. A data processing apparatus forsegmenting a small feature in a multidimensional digital array ofintensity values comprising: an input for a plurality of intensityvalues arranged along regular increments in each of a plurality ofdimensions; a memory medium for storing the plurality of intensityvalues as a multidimensional digital array; a processor configured todetect a local intensity extreme in the multidimensional digital array,to identify points along a plurality of rays originating at the totallocal intensity extreme, to identify one edge point on each ray of saidplurality of rays, said edge point associated with a maximum edge metricalong said ray, to label each point on each ray from the local intensityextreme to the edge point, and to label an unlabeled point adjacent to alabeled point if the unlabeled point has a more extreme intensity thanthe labeled point and the unlabeled point is closer than the labeledpoint to the local intensity extreme until no more unlabeled points canbe labeled; and an output for providing the labeled points forsubsequent processing.
 20. The apparatus of claim 19, wherein theplurality of intensity values arranged along regular increments in eachof a plurality of dimensions is at least one digital image, and eachpoint is a pixel.
 21. The apparatus of claim 20, wherein the digitalimage is a digitized mammogram and the small feature is amicrocalcification candidate.
 22. A computer program embodied in anon-transitory computer readable medium for performing the steps of:computing an edge metric along each ray of a plurality ofmultidimensional rays originating at a local intensity extreme,identifying a multidimensional edge point corresponding to a maximumedge metric on each said ray, labeling every point on each said ray fromsaid local intensity extreme to said edge point, and labeling anunlabeled point if the unlabeled point is adjacent to a labeled pointand the unlabeled point has a more extreme intensity than the labeledpoint and the unlabeled point is closer than the labeled point to thelocal intensity extreme.
 23. A method of labeling pixels of an image soas to designate portions of the image that are associated with anobject, the method comprising: identifying a first pixel as belonging toan object due to the first pixel having an intensity that is a localintensity extreme, wherein the first pixel is at an interior of theobject; determining that a second pixel that lies on a ray that emanatesfrom the first pixel has a maximum edge metric on the ray, wherein thesecond pixel has an intensity that is smaller in magnitude than theintensity of the first pixel; labeling the second pixel as an edge pixelthat lies on an edge of the object; determining that a third pixel thatis adjacent to the second pixel satisfies a predetermined criterionrelative to one or more of the first and second pixels; and labeling thethird pixel as belonging to the object.
 24. The method of claim 23,wherein the intensity of the first pixel is greater than the intensitiesof all pixels immediately adjacent to the first pixel.
 25. The method ofclaim 23, wherein the intensity of the first pixel is less than theintensities of all pixels immediately adjacent to the first pixel. 26.The method of claim 23, wherein the edge metric comprises a slopequotient that compares a difference between intensities of the firstpixel and a pixel that is being evaluated to a distance between thefirst pixel and the pixel that is being evaluated.
 27. The method ofclaim 23, wherein the predetermined criterion comprises the third pixelbeing disposed along a substantially straight line between the first andsecond pixels.
 28. The method of claim 23, wherein the predeterminedcriterion comprises: an intensity of the third pixel being less than anintensity of the second pixel; and a distance between the first andthird pixels being smaller than a distance between the first and secondpixels by no less than an inclusion tolerance distance.
 29. The methodof claim 23 wherein the predetermined criterion comprises an intensityof the third pixel being greater than an intensity of the second pixel.30. The method of claim 23 wherein the predetermined criterioncomprises: an intensity of the third pixel being no less than anintensity of the second pixel; and the third pixel being closer to thefirst pixel than the second pixel is to the first pixel.
 31. The methodof claim 23 wherein the predetermined criterion comprises: an intensityof the third pixel being no greater than an intensity of the secondpixel; and the third pixel being closer to the first pixel than thesecond pixel is to the first pixel.
 32. The method of claim 23 whereinthe predetermined criterion comprises: an intensity of the third pixelbeing no less than an intensity of the second pixel; and a distancebetween the first and third pixels being no more than an expansivetolerance distance greater than a distance between the first and secondpixels.
 33. The method of claim 23 wherein the predetermined criterioncomprises: an intensity of the third pixel being no greater than anintensity of the second pixel; and a distance between the first andthird pixels being no more than an expansive tolerance distance greaterthan a distance between the first and second pixels.
 34. The method ofclaim 23 wherein the predetermined criterion comprises: an intensity ofthe third pixel being no less than an intensity of the second pixel; andno less than an inclusion portion of the third pixel being on a side ofa substantially straight inclusion line closest to the first pixel, theinclusion line intersecting the second pixel and being substantiallyperpendicular to a substantially straight line that intersects the firstand second pixels.
 35. The method of claim 23 wherein the predeterminedcriterion comprises: an intensity of the third pixel being no greaterthan an intensity of the second pixel; and no less than an inclusionportion of the third pixel being on a side of a substantially straightinclusion line closest to the first pixel, the inclusion lineintersecting the second pixel and being substantially perpendicular to asubstantially straight line that intersects the first and second pixels.36. The method of claim 23, further comprising identifying as part ofthe edge of the object a fourth pixel that is immediately adjacent to atleast one pixel that is identified as part of the object and that isimmediately adjacent to at least four other pixels that are outside ofthe object.
 37. A method of labeling a subset of pixels of an image, themethod comprising: labeling pixels of an image as belonging to a firstobject that is encompassed by a first edge; labeling pixels of the imageas belonging to a second object that is encompassed by a second edge;and assembling the first and second objects into a third object that islarger than either of the first and second objects if a distance betweenthe first and second edges is no more than a join distance.
 38. Themethod of claim 37 wherein assembling the first and second objects intothe third object comprises identifying as part of the third object apixel that is disposed between the first and second edges.
 39. Themethod of claim 37 wherein assembling the first and second objects intothe third object comprises identifying as part of the third object apixel that is disposed between the first and second edges and that is nofarther than the join distance from the first edge and from the secondedge.
 40. A non-transitory computer-readable medium having instructionsstored thereon, the instructions comprising: instructions for labelingpixels of an image as belonging to a first object that is encompassed bya first edge; instructions for labeling pixels of the image as belongingto a second object that is encompassed by a second edge; andinstructions for assembling the first and second objects into a thirdobject that is larger than either of the first and second objects if adistance between the first and second edges is no more than a joindistance.
 41. A method of labeling pixels of an image so as to designateportions of the image that are associated with an object, the methodcomprising: identifying a first pixel as belonging to an object due tothe first pixel having an intensity that is a local intensity extreme,wherein the first pixel is spaced from an edge of the object;identifying as belonging to the object a second pixel that lies on afirst substantially straight line on which the first pixel also lies;identifying a third pixel as belonging to the object after havingidentified the second pixel as belonging to the object, wherein thethird pixel lies on the line at a position that is closer to the firstpixel than is the second pixel; and labeling each of the first, second,and third pixels as belonging to the object.
 42. The method of claim 41,wherein the intensity of the first pixel is greater than the intensitiesof all pixels immediately adjacent to the first pixel.
 43. The method ofclaim 41, wherein the intensity of the first pixel is less than theintensities of all pixels immediately adjacent to the first pixel. 44.The method of claim 41, wherein identifying the second pixel asbelonging to the object comprises: calculating respective slopequotients of respective differences between intensities of the firstpixel and other pixels that are intersected by the line and respectivedistances between the first pixel and the other pixels, wherein thesecond pixel is one of the other pixels; and determining that amagnitude of the slope quotient of the second pixel is larger than themagnitudes of the slope quotients for the remainder of the other pixels.45. The method of claim 41, further comprising identifying as belongingto the object at least a fourth pixel that intersects a second line thatalso intersects the first pixel but that does not intersect the secondand third pixels, the fourth pixel being identified before any pixelother than one or more of the first, second, and third pixels isidentified as belonging to the object.
 46. The method of claim 41,wherein the object comprises an edge and the second pixel forms at leasta portion of the edge.
 47. The method of claim 41, wherein the secondand third pixels are identified before any pixel other than the firstpixel is identified as belonging to the object.
 48. The method of claim41, wherein one or both of the second and third pixels have respectiveintensities that are smaller in magnitude than the intensity of thefirst pixel.
 49. The method of claim 48, wherein the intensity of thethird pixel is less than the intensity of the second pixel.
 50. Anon-transitory computer-readable medium having instructions storedthereon, the instructions comprising: instructions for identifying afirst pixel as belonging to an object due to the first pixel having anintensity that is a local intensity extreme, wherein the first pixel isspaced from an edge of the object; instructions for identifying asbelonging to the object a second pixel that lies on a firstsubstantially straight line on which the first pixel also lies;instructions for identifying a third pixel as belonging to the objectafter having identified the second pixel as belonging to the object,wherein the third pixel lies on the line at a position that is closer tothe first pixel than is the second pixel; and instructions for labelingeach of the first, second, and third pixels as belonging to the object.51. A non-transitory computer-readable medium having instructions storedthereon, the instructions comprising: instructions for identifying afirst pixel as belonging to an object due to the first pixel having anintensity that is a local intensity extreme, wherein the first pixel isat an interior of the object; instructions for determining that a secondpixel that lies on a ray that emanates from the first pixel has amaximum edge metric on the ray, wherein the second pixel has anintensity that is smaller in magnitude than the intensity of the firstpixel; instructions for labeling the second pixel as an edge pixel thatlies on an edge of the object; instructions for determining that a thirdpixel that is adjacent to the second pixel satisfies a predeterminedcriterion relative to one or more of the first and second pixels; andinstructions for labeling the third pixel as belonging to the object.